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^ Aims. We introduce and test a new and highly efficient method for treating the thermal and radiative effects influencing the energy equation in 
Q . SPH simulations of star formation. 

3*j Methods. The method uses the density, temperature and gravitational potential of each particle to estimate a mean optical depth, which 
C/j then regulates the particle's heating and cooling. The method captures - at minimal computational cost - the effects of (i) the rotational and 
vibrational degrees of freedom of H 2 ; (ii) H 2 dissociation and H° ionisation; (iii) opacity changes due to ice mantle melting, sublimation of 
dust, molecular lines, H~, bound-free and free-free processes and electron scattering; (iv) external irradiation; and (v) thermal inertia. 
Results. We use the new method to simulate the collapse of a 1 M© cloud of initially uniform density and temperature. At first, the collapse 
^ proceeds almost isothermally (T oc p 008 ; c f. Larson 2005). The cloud starts heating fast when the optical depth to the cloud centre reaches 
t^- unity (p c ~ 7 x 10~ 13 g cm -3 ). The first core forms at p c ~ 4 x 10 -9 g cm -3 and steadily increases in mass. When the temperature at the 
■ centre reaches T c ~ 2000 K, molecular hydrogen starts to dissociate and the second collapse begins, leading to the formation of the second 
(protostellar) core. The results mimic closely the detailed calculations of Masunaga & Inutsuka (2000). We also simulate (i) the collapse 
of a 1.2 M© cloud, which initially has uniform density and temperature, (ii) the collapse of a 1.2 M© rotating cloud, with an m = 2 density 
if~) perturbation and uniform initial temperature, and (iii) the smoothing of temperature fluctuations in a static, uniform density sphere. In all these 
C^*) ' tests the new algorithm reproduces the results of previous authors and/or known analytic solutions. The computational cost is comparable to 
C^- . a standard SPH simulation with a simple barotropic equation of state. The method is easy to implement, can be applied to both particle- and 
grid-based codes, and handles optical depths < r ^ 10 11 . 
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1, Introduction tend to be computationally expensive. Indeed, even the treat- 
ment of full 3D radiative transfer on a single snapshot during 
Smoothed Particle Hydrodynamics (SPH) (Lucy 1977; the evolution of a simulation is computationally quite expen- 
Gingold & Monaghan 1977) is a Lagrangian method which in- siye (Stamatellos & Whitworth 20 05; Stamatellos et al. 2005). 
vokes a large ensemble of particles to describe a fluid, by as- 
signing properties such as mass, m„ position, r„ and velocity, More often > in SPH simulations of star formation, it is stan- 
v,-, to each particle, i. Intensive thermodynamic variables like dard practice to use a barotropic equation of state, i.e. to put 
density and pressure (and their derivatives) are estimated using p = p ^ ( e -S- Bonnell 1994; Whitworth et al. 1995; Bate 
local averages (for reviews see Benz 1990, 1991; Monaghan 1998 )- The form of P <P) is chosen t0 mimic the thermody- 
1992 2005) namics of star forming gas, as revealed by computations of the 
Radiative transfer (RT) has only recently been included spherically symmetric collapse of a single, isolated protostar 
in SPH codes (Oxley & Woolfson 2003; Whitehouse & Bate ^ Boss & M y hiU 1992; Mas ™aga & Inutsuka 2000). 
2004, 2006; Whitehouse et al. 2005; Viau et al. 2005; Mayer This is not an ideal situation, (a) A barotropic equation of 
et al. 2007). These SPH-RT codes use different simplifying state is unable to account for the fact that the thermal history of 
assumptions in order to by-pass the full treatment of multi- a protostar depends sensitively on its environment, geometry 
frequency radiative transfer in 3 dimensions (a task which is and mass; for example, low-mass protostars remain optically 
not possible with current computing resources), but they still thin to their cooling radiation to higher densities than high- 

mass ones. Thus, the evolution of the density and temperature 

Send offprint requests to: D. Stamatellos cannot be approximated by a single barotropic equation for ev- 
e-mail: D.Stamatellos@astro.cf.ac.uk ery system. Even for the same system, the density and tern- 
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perature evolution away from the centre of the cloud does not 
follow the corresponding evolution at the centre of the cloud 
(Whitehouse & Bate 2006). (b) A barotropic equation of state 
is unable to capture thermal inertia effects (i.e. situations where 
the evolution is controlled by the thermal timescale, rather than 
the dynamical one). Such effects appear to be critical at the 
stage when fragmentation occurs (e.g. Boss et al. 2000). 

One of the ultimate goals of star formation simulations is to 
track the thermal history of star-forming gas. Strictly speaking, 
this requires a computational method which can treat properly, 
in 3 dimensions, the time-dependent radiation transport which 
controls the energy equation. However, this is computationally 
very expensive, significantly more expensive than the hydro- 
dynamics. We have therefore developed a new algorithm which 
enables us to distinguish the thermal behaviours of protostars of 
different mass, in different environments, with different metal- 
licities, and to capture thermal inertia effects, without treating 
in detail the associated radiation transport. The method uses the 
density, temperature and gravitational potential of each SPH 
particle (which are all calculated using the standard SPH for- 
malism) to estimate a characteristic optical depth for each par- 
ticle. This optical depth then regulates how each particle heats 
and cools. 

The paper is organized as follows. In Section [2 we present 
the new algorithm we have developed to treat the energy equa- 
tion within SPH, also describing the aspects of SPH that are re- 
quired for a complete picture of the new method. In Section [3] 
we outline the properties we adopt for the gas and dust in a 
star-forming cloud, i.e. the composition, energy equation, equa- 
tion of state, and opacity. In Section HI we present a simulation 
of the collapse of a 1 M molecular cloud; we describe in de- 
tail the different stages of the collapse, and compare our re- 
sults with those of Masunaga & Inutsuka (2000). In Section 
three additional tests are presented: (i) the collapse of a 1.2 M 
molecular cloud (Boss & Myhill 1992; Whitehouse & Bate 
2006), (ii) the collapse of a rotating 1.2 M© molecular cloud, 
with an m = 2 density perturbation (Boss & Bodenheimer 
1979; Whitehouse & Bate 2006), and (iii) the smoothing of 
temperature fluctuations in a static, uniform-density sphere 
(Spiegel 1957). Finally, in Section[6l we summarise the method 
and the tests performed, and discuss the applicability of the new 
algorithm to simulations of astrophysical systems. 



2. The method 

The key to the new method is to use an SPH particle's density, 
Pi, temperature, Tu and gravitational potential, xj/i, to estimate 
a mean optical depth, f for the SPH particle. This mean op- 
tical depth then regulates the SPH particle's radiative heating 
and cooling; in other words, it determines the extent to which 
the SPH particle is shielded from external radiation, and the 
extent to which the SPH particle's cooling radiation is trapped. 
(The gravitational potential is used here, purely because grav- 
ity is the only particle parameter which is already calculated by 
the SPH code but is not a local function of state. Therefore it 
should, in some very general sense, represent the larger-scale 
environment surrounding the SPH particle.) 



Specifically, each SPH particle is treated as if it were em- 
bedded in a spherically- symmetric pseudo-cloud (its personal 
pseudo-cloud). The density and temperature profiles of the 
pseudo-cloud are modelled with a polytrope of index n = 2, 
but the pseudo-cloud is not assumed to be in hydrostatic bal- 
ance. (We will show later that the choice of n is not critical.) 

The position of the SPH particle within its pseudo-cloud is 
not specified; instead we take a mass-weighted average over 
all possible positions (see Figs. 1 & 2). For any given po- 
sition of the SPH particle within the pseudo-cloud, the cen- 
tral density, p c , and scale-length, R , are chosen to reproduce 
the density and gravitational potential at the position of the 
SPH particle. Similarly, the pseudo-cloud's central tempera- 
ture, T c is chosen to match the temperature at the position of 
the SPH particle. (Because the pseudo-cloud is not necessar- 
ily in hydrostatic equilibrium, we cannot - in general - write 
T c = AnGihp c R 2 Q l(n + l)k B , where m is the mean gas-particle 
mass and k B is Boltzmann's constant.) 

The optical depth, r,-, is then calculated by integrating out 
along a radial line from the given position to the edge of the 
pseudo-cloud, i.e. through the cooler and more diffuse outer 
parts of the pseudo-cloud. In this way, r,- samples the different 
opacity regimes that are likely to surround the SPH particle. 
This accounts for the fact that, even if the opacity at the position 
of the SPH particle is low, its cooling radiation may be trapped 
by cooler more opaque material in the surroundings. 

Finally, f is obtained by taking a mass-weighted average 
over all possible positions within the pseudo-cloud. 

2.1. Basic SPH equations 

We use the dragon SPH code (Goodwin et al. 2004a, b). dragon 
uses variable smoothing lengths, hi, adjusted so that the num- 
ber of neighbours is exactly W NEIB = 50; it is important to have 
a constant number of neighbours to minimise numerical dif- 
fusion (Attwood et al. 2007). An octal tree is used to collate 
neighbour lists and calculate gravitational accelerations, which 
are kernel- softened using particle smoothing lengths. Standard 
artificial viscosity is invoked in converging regions, and multi- 
ple particle time-steps are used. 

The density at the position of SPH particle i is given by a 
sum over its neighbours, j, 



(1) 



where m } is the mass of particle j, hij = (hi + hj)/2, and W(s) 
is the dimensionless smoothing kernel. 

The equation of motion for SPH particle i is 



dsi 
dt 



m, 



Pi P 



ij ij 



hii 



d\j 
dt 



(2) 



Here, Pi and Pj are the pressures at the positions of particles 
i and j respectively; r^- = r ; - r 7 - ; r t j = |r^-| ; and W'(s) = 
dW/ds. Artificial viscosity is represented by the term 



n, 7 = 



Pa 



(3) 
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with a = 1 , p = 2 , Cij = (a + c-)/2 (where a = (Pi /pi) 
the isothermal sound speed of particle i), p t j = (p t + pj)/2 , 



is 



pseudo-cloud 



Pij = WW 



if VyTy <0, 

if vyr y >0, 



(4) 



and \ij = v,- - Yj. 

The second term on the righthand side of Eqn. (|2]) is the 
gravitational acceleration experienced by particle i, given by 



dVi 
dt 

where 
W 



= _ V ( Gm J Ti j 



r 3 . 



r*s =s 
Js'=0 



W(/)47r/ 2 J^. 



(5) 



(6) 



Similarly the gravitational potential at the position of particle i 
is given by 



where 

ps'-oo 

W**(s) = W*(s) + s I W r (j , )47rj , dj / . 

J s'-s 



(7) 



(8) 



In Eqns. © and ©, the sums are over all particles except i; 
the terms W*(s) in Eqn. © and W**(s) in Eqn. ([8]) represent 
kernel softening. In practice, the calculation of these gravita- 
tional terms is rendered more efficient by using a tree struc- 
ture to identify distant clusters of particles whose effect can be 
treated collectively with a multipole expansion; this reduces an 
N 2 process to an N£n[N] process, where N is the total number 
of SPH particles. 

The energy equation for particle i is 



dui Iyh I Pi Pj _ 



h 4 r - 
% r ij 



ha 



dui 
dt 



(9) 



Here U{ is the internal energy per unit mass. The first term on 
the righthand side represents compressional and viscous heat- 
ing. The second term on the righthand side is the net radiative 
heating rate; this paper is primarily concerned with the evalua- 
tion of this term. 

2.2. Calibrating the pseudo-cloud 

Suppose that SPH particle i is embedded at radius R = £R , 
in a pseudo-cloud with central density p c , scale-length R Q , and 
poly tropic index n (see Figs. 1 & 2). f is thus a dimensionless 
radius, and p c and R Q are chosen so as to reproduce - at this 
radius - the actual density and gravitational potential of the 
SPH particle, i.e. 



-AnG Pc R 2 m = *i 



(10) 
(11) 




possible positions of 
the SPH particle 
inside its pseudo-cloud 
(dashed-circles) 



Fig. 1. Schematic representation of the pseudo-cloud around an 
SPH particle.The location of the SPH particle inside its pseudo- 
cloud is not specified. 

Here #(f) is the Lane-Emden Function for index n 
(Chandrasekhar 1939), 



d6 

m = ~^d~^ + m > 



(12) 



and f B is the dimensionless boundary of the polytrope (i.e. the 
argument of the smallest zero of 0(f)) 

If we fix n (and hence the forms of #(f ) and 0(f)), and we 
pick an arbitrary value for f (modulo that it must be within the 
pseudo-cloud, i.e. f < f B ), then we obtain 



R = 



4nG Pi </>({) 



1/2 



(13) 
(14) 




Fig. 2. Density and temperature profiles for a polytropic 
pseudo-cloud with n = 2. The SPH particle could be located 
anywhere in the cloud (solid and dashed line circles). 
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In an analogous manner we chose the central temperature 
of the pseudo-cloud so as to reproduce - at radius R = £R - 
the actual temperature of the SPH particle (see Fig. 2), 



T c = TiGT l (g). 



(15) 
(16) 



The column-density on a radial line from this radius to the 
boundary of the pseudo-cloud is then given by 



r 



1/2 



r 



(17) 



To obtain the pseudo-mean column-density, we take a mass- 
weighted average of £,•(£) over all possible dimensionless radii, 
f , i-e. 



2; 



~ in [4nG\ 



(?b) 



-1 



r 



(18) 



where [- ^ ^(f B )] I s me tota l dimensionless mass of the poly- 
trope, 0"(£)£ 2 d£ is the dimensionless mass element between £ 
and £ + and 



,d0 



4 f^ B ,=f 

J^=0 Jf'=f 



0© 



1/2 



As an indication of how insensitive the results are to the 
choice of n, we note that £ = 0.376, £ 15 = 0.372, £ 2 = 
0.368, 4 5 = 0.364, and £ = 0.360 . Since for protostars 
which are close to equilibrium - for example those undergoing 
quasistatic (i.e. Kelvin-Helmholtz) contraction - the poly tropic 
exponent is likely to fall in the range 4/3 to 5/3, we adopt 
n = 2 (corresponding to a poly tropic exponent of 3/2). 

We calculate the pseudo-mean optical depth in the same 
way as the pseudo-mean column-density. If the Rosseland- 
mean opacity is a function of density and temperature only, 
K R (p, T), the radial optical depth from radius R = £R Q to the 
boundary of the pseudo-cloud is 



(PcW.W)) p^(g)R dg 



4 0(f) 



1/2 



0(f) 



0(O 



0(0) 



0(O 



0(C) 



0(O 



dg, (20) 



and the mass-weighted pseudo-mean optical depth is 

i-i . ... - -.1/2 r?=t B 



Ti = 



d^ 
0(g) 



K R\Pi 



0(O 



4nG 

0(g)] 



x 



0(O 



e n (g)dg 



no 



1/2 



At first sight it might appear that the double integral in Eqn. 
(I2T1) will have to be evaluated on-the-fly for every SPH particle, 
at every time-step. In fact, we can define a pseudo-mean mass 
opacity 



(22) 



which - once n has been fixed - is simply a function of pi and 
Ti. It can therefore be evaluated in advance, once and for all 
time, and stored in a dense look-up table for subsequent refer- 
ence and interpolation. For the point (p, T) in the table, 



Kip,T) = 




dgfdg (23) 



The physical interpretation of this pseudo-mean opacity, 
K R (p, T), is fundamental to the method. Although formally 
k R (p, T) only depends on the local density and temperature, 
when it is multiplied by the pseudo-mean column-density, £/, 
it gives a pseudo-mean optical depth, f which allows for the 
fact that radiation absorbed or emitted by particle / has to pass 
through surrounding material which will in general have differ- 
ent density and temperature, and hence different opacity. For 
example, an SPH particle whose local Rosseland-mean opacity, 
K R (p, T), is low because its density and temperature fall in the 
opacity gap, will have a larger pseudo-mean opacity, k R (p, T)\ 
this simply reflects the fact that this SPH particle may still be 
well insulated by cooler material in its surroundings which has 
much higher opacity because it contains dust. 

2.3. Radiative heating and cooling 

The net radiative heating for SPH particle i is given by 
4cT ca (T 4 o (n)-Tf) 



dui 
dt 



^K R (puTi) + K^(p h Ti) 



(24) 



where <x SB is the Stefan-Boltzmann constant, k R (p,T) is the 
pseudo-mean opacity defined in Section [2T2l and K p (p, T) is the 
Planck-mean opacity. 

The positive term in Eqn. (l24l) - the one involving T^fa) 
- represents radiative heating due to the background radiation 
field with effective temperature r o (r/). This term ensures that 
the SPH particle does not cool radiatively below r o (r ; ). In a 
simulation which includes stars - either pre-existing, or formed 
as an outcome of the simulation; and with luminosities L+ and 
positions r* - we set 

U 



7>) = (10 K) 4 + J] 



16 ti cr |r-r*| : 



(25) 



(21) 



The negative term in Eqn. (l24l) - the one involving Tf - rep- 
resents radiative cooling of SPH particle i. If Tf » T^(ri), we 
can neglect the heating term and consider two limiting regimes: 

(i) If £^ K R (pi, Ti) <?c K~ l (pu Ti), we are in the optically thin 
cooling regime and Eqn. (l24t approximates to 

dw 



-4cr SB T t K p (p h Ti), 



(26) 
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in exact agreement with the definition of the Planck-mean 
opacity. 

(ii) If T?k R (pi,Ti) » K~ l (pi,Ti), we are in the optically 
thick cooling regime and Eqn. (l24l) approximates to 



dui 
dt 



Ca SB T j 



(27) 



where c is the speed of light, a SB is the radiant energy density 
constant, and we have obtained the second expression by sub- 
stituting 4<r SB = ca SB and % \(pu Ti) = f;. 

To see that this is just the diffusion approximation, suppose 
that the pseudo-cloud has pseudo-mass M t and pseudo-radius 
Ri ~ (3M i -/47r£0 1/2 • Eqn. ([27]) then reduces to 



dui 
dt 



rad,/ 



(28) 



Mi t dm 

where £/ ra d,z ~ 4nR^a SB Tf /3 is the total radiant energy in the 
pseudo-cloud and t m ^,i ~ RiTi/c is the timescale on which radia- 
tion diffuses out of the pseudo-cloud (cf . Masunaga & Inutsuka 
1999; Whitworth & Stamatellos 2006). 



2.5. Method implementation 

The method is easy to implement. In practice we do the follow- 
ing, for each SPH particle i, at each timestep: 

1. Calculate the pseudo-mean column-density % from the 
density, pi and gravitational potential, using Eqn. (fT8l) . 
(In a simulation which includes stars, we must neglect their 
contribution to fa.) 

2. Calculate the pseudo-mean opacity, k R (pi,Ti) and the 
Planck-mean opacity, /c p (p/, Ti), by interpolation on a look- 
up tableQ. 

3. Calculate the compressive plus viscous heating rate, 
dui/dt\ RYDRQ , using Eqn. (|29k and the radiative heating rate, 
dui/dt\ RAD , using Eqn. (l24l) . 

4. Calculate the equilibrium temperature, T Qq j, using 
Eqn. (l30l) and the thermalization timescale, ^hemv, using 
Eqn. (ED. 

5. Update the internal energy, U[, using Eqn. (l32l) ; and hence 
also advance the temperature, TV 



2.4. Quasi-implicit scheme 

In order to avoid very short time-steps, we use the following 
scheme to update the internal energy, ui. From SPH we know 
the net compressive plus viscous heating rate, 



dui 
dt 



1^ 



(29) 



We can therefore calculate (a) the equilibrium temperature T Qq j 
for each particle from 



dui 
dt 



4cr SB [7>)-7^J 
2^k R (pi,T QqJ ) + K- l (pi,T^i) 



0; 



(30) 



(b) the equilibrium internal energy, u eq j = uipi, T Qq ^)\ and (c) 
the thermalization timescale, 



f i ) J duj 

nherm,i — j^eq,/ Un \ 



dt 



dUi 

dt 



We then advance Ui through a time step At by putting 



Ui{t+At) = Ui(t) exp 



■At 



^therm,/ 



+ Meq,i S 1 - exp 



-A? 



^the 



(31) 



• (32) 



If At <?c ^therm,/, we are in a situation where thermal inertia 
effects are important, and Eqn. (l32l) approximates to 



+ AO ^ n f (r) + [w eq ,*- - 



^therm,;' 



(33) 



i.e. in one dynamical timestep, Ar, the gas can only relax a little 
towards thermal equilibrium. 

On the other hand, if At » *thenn,i» Eqn. (l32l) approximates 

to 



+ AO ^ w 



eq,« 5 



(34) 



i.e. thermal processes are occurring much faster than dynamical 
ones, and the gas is always close to thermal equilibrium. Using 
the above procedure we capture thermal inertia effects, whilst 
avoiding the use of very small timesteps. 



2.6. Limitations of the method 

Although the method is very efficient, it evidently has limita- 
tions, in particular: 

1. Because the diffusion approximation is applied here glob- 
ally (to the whole pseudo-cloud), the method cannot cap- 
ture in detail the local nature of radiative heating and cool- 
ing in the optically thick regime (i.e. the fact that in reality 
fluid elements exchange heat directly with other fluid ele- 
ments, within a few photon mean-free-paths). 

2. Because the pseudo-cloud Ansatz predicates a spherical 
polytropic cloud, the method works best for configurations 
which approximate to spherical symmetry. For example in 
simulations of disc fragmentation it handles the conden- 
sations better than the background disc. Notwithstanding 
this, even in an unperturbed disc the method is reasonably 
accurate, as shown by its performance of the Hubeny test 
(Stamatellos & Whitworth, in preparation). 



3. Gas and dust properties 

3.1. Gas-phase chemical abundances 

Although metals make essential contributions to the opacity, 
they make very little contribution to the equation of state. 
Therefore, for the purpose of treating the gas-phase chemistry, 
we assume that the gas is 70% hydrogen and 30% helium by 
mass: X = 0.7,7 = 0.3,Z = 0.At low temperatures, hydrogen 
is molecular, but as the temperature increases it becomes disso- 
ciated and then ionised. At low temperatures, helium is neutral 
atomic, but as the temperature increass it becomes ionised, first 
to He + , and then to He ++ . The relative abundances of these con- 
stituents depend on the density, pu and the temperature, Ti, and 
are calculated using Saha equations (e.g. Black & Bodenheimer 



1 Tabulated pseudo-mean opacities and internal energies (see next 
Section) can be obtained by contacting D.Stamatellos@astro.cf.ac.uk 
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1975), with the simplifying assumption that the dissociation of 
H 2 is complete before ionization of H° begins; and similarly, 
that the ionization of He is complete before the ionization of 
He + begins. 




Temperature (K) 



Fig. 3. The variation of the mean molecular weight with den- 
sity and temperature. Isopycnic curves are plotted from p = 
10" 18 gem -3 to p = 1 gem -3 , every two orders of magnitude 
(bottom to top). 




Temperature (K) 



Fig. 4. The variation of the specific internal energy with den- 
sity and temperature. Isopycnic curves are plotted from p = 
10~ 18 gcm~ 3 to p = 1 gem -3 , every two orders of magnitude 
(from top to bottom). The rotational degrees of freedom of H2 
are excited around 100 to 200 K, H2 starts to be dissociated 
around 1, 000 to 10, 000 K, and H° starts to be ionised around 
4,000 to 40, 000 K. 



where 



3.2. The equation of state 

If we define y = nu° /2nu 2 to be the degree of dissociation of 
hydrogen, x = nu+ Mh° to be the degree of ionization of hy- 
drogen, z\ = ^He + / n n&° to be the degree of single ionisation of 
helium, and zi - ^He ++ / n n& + to be the degree of double ionisa- 
tion of helium, then the mean molecular weight is given by 



fii = ji{pu Ti) 



X Fl -1 
(1 + y + 2xy) - + (1 + zi + Z\Zi) - . 



Note that (x,y,zi,Z2) must be evaluated afresh for each SPH 
particle; the index i has been dropped purely for simplicity. The 
variation of the mean molecular weight with density and tem- 
perature is shown in Fig. [51 

For densities up to ~ 0.03 g cm -3 the ideal gas approxima- 
tion holds, and hence the gas pressure is 



Pi = 



pi K T i 

{limn 



(36) 



un 2 



^He 



^H 2 DISS 



X(l-y) 



uu = Xy (1 + x) 



3KT t 



2m h 



2m H 



Y(\+z\ +Z1Z2) 



3k B Tj 
8mR 



Xy 



H 2 Diss 



(35) u m0N = Xxy 



2m H 

^HION 



^He ION 



^He + ION - 



Yzi(l-z 2 ) 



Jh, 



elON 



4m H 

-^He+ION 



4m H 



(38) 
(39) 
(40) 
(41) 
(42) 
(43) 
(44) 



Here, Dh 2 diss = 4.5 eV is the dissociation energy of H2; 
Jhion = 13.6 eV, J He iON = 24.6 eV and J He +iON = 54.4 eV 
are the ionisation energies of H°, He and He + , respectively; 
and the function 



3.3. Specific internal energy of the gas 

The specific internal energy (energy per unit mass) of an 
SPH particle i is the sum of contributions from molecular, 
atomic and ionised hydrogen, atomic, singly ionised and dou- 
bly ionised helium, and the associated dissociation and ionisa- 
tion energies, 

Ui = Uu 2 + Uu + Uuq + Wh 2 DISS + ^HION + ^HelON + ^He+ION, (37) 



, r y /^rot\ rsrj, x , / ^vib \ ex p(^vm/^) 



(45) 



with r ROT = 85.4 K and T wlB - 6100 K accounts for the rota- 
tional and vibrational degrees of freedom of H 2 . The function 
f(Ji) depends on the relative abundances of ortho- and para- 
H 2 ; we assume a fixed ortho-to-para ratio of 3:1. The variation 
of the specific internal energy with density and temperature is 
shown in Fig.lU 
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3.4. Opacity 

In the present work we do not distinguish between the 
Rosseland-mean and Planck-mean opacities; we use the 
parametrisation proposed by Bell & Lin (1994) for both, i.e. 

K R (p,T) = * p (p,D = K p a T b . (46) 

Here (K ,a,b) are constants which depend on the dominant 
physical process contributing to the opacity in different regimes 
of density and temperature (see Table Q] and Fig. [5]). 

The opacity at low temperatures is dominated by icy dust 
grains. At T ~ 150 K the ices evaporate and the opacity is due 
to metal grains up to T ~ 1 , 000 K, when the metal grains start 
to evaporate. The opacity drops considerably in the tempera- 
ture range from T ~ 1,000 K to T ~ 2,000 K, as it is too 
hot for dust to exist, and too cool for H" to contribute, so the 
opacity is mainly due to molecules; this region of low opacity 
is sometimes referred to as the opacity gap. The opacity starts 
to increase again above T ~ 2, 000 K due to H" absorption and 
then decreases again above T ~ 10 4 K, when free-free transi- 
tions take over. At very high temperatures, electron scattering 
delivers an approximately constant opacity. 

At low temperatures, T < 2,000K, the Bell & Lin 
parametrisation agrees well with the Rosseland-mean dust 
opacity calculated by Preibisch et al. (1993). Similarly, at high 
temperatures, T > 2, 000 K, it agrees well with the Rosseland- 
mean gas opacities calculated by Alexander & Ferguson (1994) 
and Iglesias & Rogers (1996). 

Eqn. (l46t gives local Rosseland- and Planck-mean opaci- 
ties. To calculate the pseudo-mean opacity used in Eqn. (|24k 
we have to convolve this opacity with polytropic density and 
temperature profiles according to Eqn. (l23t . In Fig. [6] we 
present the pseudo-mean opacity computed in this way, using a 
polytropic index n = 2. We reiterate that the choice of n affects 
the computed pseudo-mean opacities only weakly. 



Table 1. Opacity law parameters (from Bell & Lin 1994) 





Dominant opacity component 
or physical process 


K 


a 


b 


1 


Ice grains 


2 x 10- 4 





2 


2 


Evaporation of ice grains 


2x 10 16 





-7 


3 


Metal grains 


0.1 





1/2 


4 


Evaporation of metal grains 


2 x 10 81 


1 


-24 


5 


Molecules 


io- 8 


2/3 


3 


6 


H~ absorption 


10 -36 


1/3 


10 


7 


bf and ff transitions 


1.5 x 10 20 


1 


-5/2 


8 


Electron scattering 


0.348 









4. The collapse of a 1-Mq molecular cloud 

The first test of our new method for treating the energy equation 
in SPH is to simulate the collapse of a 1 M Q molecular cloud, 
which initially is spherical with radius R = 10 4 AU and uni- 
form density po = 1.41 x 10" 19 g cm -3 . We set the background 
radiation temperature to T Q (r) = 5 K This problem has been 




10 10 2 10 3 10 4 10 5 10 6 

Temperature (K) 

Fig. 5. The variation of the local Rosseland-mean opacity with 
density and temperature. Isopycnic curves are plotted from 
p = 10" 18 gcm -3 to p = 1 gem -3 , every two orders of mag- 
nitude (from bottom to top). The opacity gap is evident at tem- 
peratures ~ 1, 000 to 3, 000 K, over a wide range of densities. 




Temperature (K) 

Fig. 6. The variation with density and temperature of the 
pseudo-mean opacity. Isopycnic curves are plotted as in Fig. 
For comparison the local opacity at density p = 10" 6 g cm -3 is 
also plotted (dashed line). 



investigated by Masunaga & Inutsuka (2000) using a code that 
treats the hydrodynamics in 1 dimension (i.e. assuming spher- 
ical symmetry) and the radiative transfer exactly in 3 dimen- 
sions (i.e. by solving the angle-dependent and frequency de- 
pendent radiation transfer equation). It therefore constitutes a 
stiff test for our new method to reproduce their results. For the 
simulation presented here we use 2 x 10 5 SPH particles. 

4.1. Cloud collapse and the formation of the first and 
second cores 

The evolution of the cloud is followed up to density p ~ 
10" 3 gem -3 and temperature T ~ 10 4 K, i.e. a difference from 
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(8) Electron scattering 



(7) B-F/F-F transitions 




(1) Ice grains 



I i i i i i i i i i i i i i i i I 

10 -18 1Q -17 10 -16 1Q -15 1Q-14 10 -13 lO"^ 10" 11 10" 10 10" 9 10" 8 10" 7 10" 6 10" 5 1 0" 4 10" 3 10-2 

Density (g cm -3 ) 

Fig. 7. Evolution of the central density and central temperature of the collapsing cloud (thick red line). For comparison, the 
dashed black line shows the Masunaga & Inutsuka (2000) simulation, and the dotted black lines delineate the different opacity 
regimes (see Fig. [5]). The thin solid lines are the loci for different degrees of H2 dissociation (bottom set of blue lines, y = 
0.01, 0.2, 0.4, 0.6, and 0.8 , respectively) and of H° ionisation (top green line, x - 0.01). The results of our model are very close 
to the simulation of Masunaga & Inutsuka (2000). Differences at high densities are attributable to our using different opacities. 



10 




10 _i6 : — , , , . 1 . , , , 1 , , . , l_ 

1.02 1.03 1.04 1.05 

time (t°) 



Fig. 8. Evolution of the central density with time, given in units of the initial free fall time % = 1.781 x 10 5 yr. The times of the 
formation of the fist and second core are also marked. The red squares correspond to the Masunaga & Inutsuka (2000) simulation. 
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the initial conditions of about 16 orders of magnitude in den- 
sity, and 3 orders of magnitude in temperature. 

As long as the density is below ~ 10~ 12 gcm~ 3 , the 
temperature increases slowly with increasing density. For 
10" 18 g cm -3 ^ p ^ 10" 13 g cm -3 we can approximate this with 



5K 



0.08 



10~ 18 gcnr 



(47) 



(cf. Larson 1973; Low & Lynden-Bell 1976; Masunaga & 
Inutsuka 2000; Larson 2005). 

The cloud core starts heating more rapidly when it becomes 
optically thick. This continues until the the temperature reaches 
T ~ 100 K at density p ~ 3 x 10" 11 g cm -3 , when the rotational 
degrees of freedom of H2 start to get excited, and hence the 
temperature increases more slowly with density (see the kink 
around T ~ 100 K in Fig. [7]). As the temperature increases 
further, the thermal pressure starts to decelerate the contrac- 
tion, and the first core is formed (Larson 1969; Masunaga et al. 
1998; Masunaga & Inutsuka 2000; Whitehouse & Bate 2006) 
at t = 1.048 tff , where % is the free-fall time at the start of 
collapse, i.e. t So = [3n/(32 Gp )] 1/2 = 1.781 x 10 5 yr (see 
Fig.©. 

The first core grows in mass, contracts, and heats up un- 
til the temperature reaches T ~ 2, 000 K, when H2 starts to 
dissociate. Consequently the compressional energy delivered 
by contraction does not all go to heat the core; instead some 
of it goes into dissociating H2, and the second collapse starts. 
This second collapse proceeds until almost all of the molecular 
hydrogen at the centre has been dissociated. When the den- 
sity reaches p ~ 10" 3 gem -3 and the temperature rises above 
T ~ 10, 000 K, the collapse again decelerates and the second 
core (i.e. the protostar) is formed (Larson 1969; Masunaga & 
Inutsuka 2000) at t = 1.052 t$ . At first, the second core pul- 
sates (see Fig[8]and Larson 1969), but eventually it settles down 
into quasistatic contraction. Due to computational constraints, 
the evolution is not followed further. 

The evolution of the core in our simulation is very similar 
to that obtained by Masunaga & Inutsuka (2000), as shown in 
Fig El Differences at densities > 5 x 10" 6 g cm" 3 are attributable 
to the different opacities we use. The timescales in our simu- 
lation are also very similar to those obtained by Masunaga & 
Inutsuka (2000), as is shown in Fig. [U where the evolution of 
the density at the centre of the cloud is plotted against time. 
The times computed from our simulation fit well with the times 
from the Masunaga & Inutsuka (2000) simulation if we syn- 
chronise the two simulations at density p = 4.34xl0" 13 g cm -3 , 
to avoid discrepancies due to small differences in the initial 
conditions at the onset of the collapse. 

4.2. Snapshots during the cloud collapse 

In order to describe the evolution away from the centre of the 
cloud, and to make a more detailed comparison with the re- 
sults of the Masunaga & Inutsuka (2000) simulation, we fo- 
cus on eight representative instants during the cloud evolution. 
Critical parameters at these instants are listed in Table [2] The 
central densities have been chosen so as to match those used by 



Table 2. Colour code, elapsed time (t, measured from the be- 
ginning of the simulation and in units of the initial freefall 
time), central density (p c ), and central temperature (T c ), for 
the instants illustrated in Figs.l9ltofT2l 



Label 


Colour 


t/tffo 


p/gem 3 


T/K 


1 




5.22 x 10" 2 


3.16 x 10" 19 


5.0 


2 


black 


8.43 x 10" 1 


1.68 x 10" 18 


5.7 


3 


green 


1.03038936 


5.93 x 10" 16 


8.0 


4 


cyan 


1.04755301 


1.09 x 10" 10 


140 


5 


magenta 


1.04940519 


1.59 x IO- 9 


440 


6 


black 


1.05153974 


1.00 x 10" 7 


1270 


7 


green 


1.05155811 


3.30 x 10- 4 


5530 


8 


red 


1.05155816 


2.31 x 10" 3 


10210 



Masunaga & Inutsuka (2000) for the same purpose (see their 
Fig. 1 and their Table 1). 

Fig. [9] shows the run of temperature against density at the 
instants defined in Table [2] During the early stages of the col- 
lapse the variation of temperature with density mimics the evo- 
lution of the central temperature and density (see the green 
points in Fig.|9j points representing previous instants are over- 
lapped). At later stages, the regions around the centre start heat- 
ing at lower densities, as in the simulation of Whitehouse & 
Bate (2006). 

In Figs. [TOl through fT2l we present the density, temperature 
and radial infall velocity profiles at different instants during the 
evolution of the cloud. These profiles are very similar to those 
reported by Masunaga & Inutsuka (2000). The velocity profiles 
(Fig. [12]) clearly show the formation of an accretion shock at 



CD 

Z 10 3 

■+■> 

CO 

CD 



(8) Electron scattering 



(7) B-F/F-F transitions 



(6) H" absorption 



(5) Molecules 
(4) Metal grain evaporation 




■(2) Ice evaporation 



io-iBio-ino-'oio-^o-^io-^o-^o-^io-^io- 8 io- e 10- 7 io-« 10- 5 10-* io- 3 o.o] 
Density (g cm -3 ) 

Fig. 9. The run of temperature against density at different in- 
stants during the cloud evolution (instants 2 to 8 in Table [2]). 
The region around the centre of the cloud heats at lower densi- 
ties than the centre of the cloud. The density increases as time 
evolves (black, green, cyan, magenta, black, green, red). For 
reference, we also plot the evolution of the temperature and 
density at the centre of the cloud in our simulation (lower solid 
black line) and in the Masunaga & Inutsuka (2000) simulation 
(dashed black line). 
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(AU) 



Fig. 10. Density profiles at different instants during the cloud 
evolution (instants 2 to 8 in Table [2]). The colour coding is the 



same as in Fig. [9] 




Fig. 11. Temperature profiles at different instants during the 
cloud evolution (instants 2 to 8 in Table [2]). The colour cod- 
ing is the same as in Fig.O 



the boundary of the first core at radius R ~ 3 to 5 AU. There 
is also an accretion shock at the boundary of the second core, 
initially at a radius of R ~ 0.003 AU, but later expanding to 
R ~ 0.01 AU (cf. Larson 1969). 

4.3. Convergence 

We repeat this simulation using different numbers of SPH par- 
ticles, TV, to check for convergence. We use TV ^ 2 x 10 4 , 5x 
10 4 , 10 5 , and 2 x 10 5 (Fig. [13]). The run of central temperature 
against central density is almost identical for different numbers 
of particles, and the results are fully converged up to densities 
p ~ 0.003 g cm -3 with N ^ 10 5 . Currently available supercom- 
puting facilities allow SPH simulations of star formation with 
up to 3 x 10 7 particles, and so the convergence condition quoted 
above can easily be met. 



(8) Electron scattering 



(5) Molecules 

(4) Metal grain evaporation 



(7) B-F/F-F transitions 



(6) H" absorptio 



(3) Metal grains 




(2) Ice evaporation 



(1) Ice grains 



10 -ie 10 -n 10 -is lo-is 10 -i 



10 -i3 10 -is 10 -n lO-m 10 -8 10-8 iQ-t iq-b 10 -s 

Density {g cm" 3 ) 



10-* lO- 3 lO- 8 



Fig. 13. The run of central temperature against central den- 
sity during the evolution of a 1 M molecular cloud, simu- 
lated using different numbers of SPH particles: N - 2 x 
10 4 (green), 5 x 10 4 (blue), 10 5 (cyan), and 2 x 10 5 (red). 
The results converge for N t 10 5 . The black dashed line repre- 
sents the Masunaga & Inutsuka (2000) simulation. 




' i i i i i i i ; 

lO" 3 0.01 0.1 1 10 10 a 10 3 10* 

r (AU) 

Fig. 12. Radial infall velocity profiles at different instants dur- 
ing the cloud evolution (instants 2 to 8 in Table [2]). The colour 
coding is the same as in Fig. [9] 



5. Additional tests 

5.1. Boss &Myhill (1992) 

The second test of our new method is to simulate the evo- 
lution of a 1.2 M© cloud with uniform initial density p = 
1.7 x 10~ 19 gcm~ 3 , uniform initial temperature T = 10 K, and 
initial radius R = 1.5 x 10 17 cm, as originally investigated by 
Boss & Myhill (1992). This problem has recently been revis- 
ited by Whitehouse & Bate (2006), using SPH with flux-limited 
diffusion. In our simulation we have N — 1.5 x 10 5 SPH par- 
ticles. The evolution of the cloud is very similar to the evolu- 
tion already described in Section |U Fig. [14] compares the run 
of central temperature against central density which we obtain, 
with that obtained by Whitehouse & Bate (2006). There are 
three small differences, (i) In our simulation, the cloud starts 
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i i i i i i i i i i i i i i i 

10-' 8 10-' 7 10-'« 10-' 6 lO" 14 lO" 13 10-' 8 10" 11 10-'° IO- 9 lO- 8 10- 7 10-« 10- 6 10~ 4 lO" 3 10" 8 

Density (g cm -3 ) 

Fig. 14. Evolution of the central density and central tempera- 
ture for the Boss & Myhill (1992) test problem. The dashed 
line corresponds to the Whitehouse & Bate (2006) simulation, 
and the dotted lines define the different opacity regimes (see 
Fig. [5]). The results of our model are very close to the simu- 
lation of Whitehouse & Bate (2006). Differences at densities 
> 5 x 10" 6 g cm -3 are attributable to the use of different opaci- 
ties. 

heating before it becomes optically thick (as reported also by 
Masunaga & Inutsuka (2000) in a similar test). In contrast, the 
Whitehouse & Bate (2006) cloud remains strictly isothermal 
during this phase, (ii) In the Whitehouse & Bate (2006) sim- 
ulation the centre of the cloud heats up more rapidly at den- 
sities p ~ 10" 11 to 10~ 6 gcm~ 3 , i.e. at lower densities than 
in our simulation. However, Masunaga & Inutsuka (2000) in 
a similar test also find lower temperatures than Whitehouse & 
Bate (2006) in this regime, (iii) There are differences at densi- 
ties p ^ 5 x 10~ 6 gcm~ 3 , which are attributable to the use of 
different opacities. These differences are small and the overall 
evolution of the cloud is similar in the two simulations. 



I I I I I I I I I I I I I I I i i linn 

10 -18 10 -17 10 -16 10 -16 1Q-14 1Q -13 10-12 10-H 10-10 lQ-9 1Q -8 lQ -7 IQ-6 1Q -6 1Q -4 JQ-3 IQ-2 

Density (g cm -3 ) 

Fig. 15. Evolution of the density and temperature at the densest 
part of a collapsing, rotating molecular cloud. The dashed line 
corresponds to the Whitehouse & Bate (2006) simulation of the 
same problem. The results of our simulation are very close to 
those of Whitehouse & Bate (2006). Differences at densities 
p > 5 x 10" 6 g cm -3 are attributable to different opacities. 



the xy-plane at three instants during the evolution. The com- 
ponents of the binary system are connected by a bar which 
subsequently fragments. If the collapse were isothermal, this 
bar should not show any tendency to fragment (e.g. Truelove 
et al. 1998; Klein et al. 1999; Kitsionas & Whitworth 2002). 
However, Bate & Burkert (1997) show that if the gas is allowed 
to heat, but the heating happens at sufficiently high densities, 
p ^ 0.3 x 10" 13 g cm -3 , then the bar does fragment. In our sim- 
ulation, the gas in the bar starts to heat up rapidly only when 
the density reaches p ~ 7 x 10" 13 gem -3 . Therefore the frag- 
mentation of the bar is consistent with the predictions of Bate 
& Burkert (1997). However, we note that Whitehouse & Bate 
(2006) do not report any bar fragmentation. 



5.2. Boss & Bodenheimer (1979) 



5.3. Thermal relaxation 



The third test of our new method is to simulate the evolu- 
tion of a rotating 1.2 M© cloud; the cloud initially rotates as 
a solid body with angular velocity Q = 1.6 x 10" 12 rad s" 1 , and 
so the ratio of rotational-to-gravitational energy is f3 = 0.26; 
the density includes an m = 2 perturbation, i.e. p = 1.44 x 
10" 17 gcm" 3 [l + 0.5 cos(20)], where is the azimuthal angle 
in the plane perpendicular to the axis of rotation; the initial ra- 
dius of the cloud is R = 3.2 x 10 16 cm and the temperature 
is initially uniform at T = 12 K; this problem was originally 
investigated by Boss & Bodenheimer (1979), and has recently 
been revisited by Whitehouse & Bate (2006). In our simulation 
we have N - 1.5 x 10 5 SPH particles. Fig. [T5l compares the 
evolution of the density and temperature of the densest part of 
the cloud in our simulation, with that obtained by Whitehouse 
& Bate (2006). The differences are again only small. 

The result of the collapse is a binary with separation S ~ 
500 AU. In Fig. [16] we plot the density and the temperature on 



Finally, we test the time-dependence of our new method by 
simulating the relaxation of temperature fluctuations in a static 
sphere with uniform density p = 10" 19 g cm -3 and radius R = 
10, 000 AU. We assume an equilibrium temperature of T = 
10 K and an initial temperature perturbation of the form AT = 
AT Q sin(kr)/kr (Masunaga et al. 1998; Spiegel 1957), where 
AT Q = 0.15 K is the amplitude of the perturbation and k = 
7r/(2500 AU) is its characteristic wavenumber. Masunaga et al. 
(1998) have shown that at subsequent times the temperature 
should be 



T(r,t) = T + Ar ^^-"» 



kr 



where 



oj(k) = y 



(48) 



(49) 
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Fig. 16. Three instants from our simulation of the Boss & Bodenheimer test problem, at t\ - 0.02315 Myr, ti - 0.02343Myr, and 
t3 - 0.025 Myr (left to right). We plot the logarithmic density (top) and the logarithmic temperature (bottom), on the xy-plane, 
i.e. the plane perpendicular to the rotation axis. The central densities and central temperatures of the southern condensation are 
p ls = 1.2 x 10" 12 gcm- 3 , r is = 25K,p 2S = 1.6 x 10- 10 gcnr 3 , T 2S = 150K,p 3S = 3.8 x 10- 3 gcnr 3 , T 3S = 10800 K. The 
corresponding central densities and central temperatures of the northern condensation are p 1N = 1.7 x 10" 13 g cm -3 , T lN = 18 K, 
p 2N =2.4xlO- 12 gcm- 3 ,r 2N =35K,p 3N =6.3xl0- 10 gcm- 3 ,r 3N = 780K. 



is the relaxation rate, 



k q is the opacity at the equilibrium temperature, and cy is the 
heat capacity of the material. 

In Figs. [T7] to [20] we present the radial temperature pro- 
files at different instants during the relaxation towards equi- 
librium for spheres having different optical depths, r = 
0.1, 1, 10, and 100. The simulation results approximate well 
to Eqn. (HU). 

The relaxation rates are also reproduced well. In Fig. [21] 
we present the dispersion relation, i.e. the relaxation rate a>(k) 
for different values of the ratio K Q /k. We plot the relaxation 
rates for all the SPH particles that represent the uniform den- 
sity sphere, calculated at seven different snapshots during the 
temperature relaxation (dots that saturate to form lines). The 



theoretical values (calculated from Eq.[49]) are also plotted (red 
line and squares). 

6. Summary 

We have developed a new method to treat the influence of ra- 
diative transfer on the energy equation in SPH simulations of 
star formation. The method uses the density, temperature and 
gravitational potential of each particle to make an educated es- 
timate of the mean optical depth which regulates its heating and 
cooling. It can treat both the optically thin and optically thick 
regimes. 

The energy equation takes account of heating by compres- 
sion or cooling by expansion (i.e. PdV work); viscous dis- 
sipation; external irradiation; and radiative cooling. In situ- 
ations where the thermal timescale is much longer than the 
dynamical timescale, the resulting thermal inertia effects are 
captured properly. Conversely, where the thermal timescale is 
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Fig. 17. Thermal relaxation of a static, uniform sphere of opti- Fig. 19. Same as Fig. [T7]but for a sphere of optical depth r 
cal depth r = 0.1; seven instants are plotted, showing the tern- 10. 
perature relaxing to its equilibrium value T eq = 10 K. 
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Fig. 18. Same as Fig. [171 but for a sphere of optical depth r = 1 . 
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Fig. 20. Same as Fig. [T7]but for a sphere of optical depth r 
100. 



much shorter than the dynamical timescale, we avoid very short 
timesteps, essentially by assuming thermal equilibrium. 

The equation of state and the internal energy take account 
of (i) the rotational and vibrational degrees of freedom of H2, 
and (ii) the different chemical states of hydrogen and helium 
(cf. Black & Bodenheimer 1975; Boley et al. 2007) 

A simple parametrisation of the frequency averaged opac- 
ity is used (Bell & Lin 1994), which reproduces the basic fea- 
tures of more sophisticated opacity models (e.g. Preibisch et al. 
1993; Alexander & Ferguson 1994). This parametrisation ac- 
counts for the effects of dust sublimation, molecules, H", free- 
free transitions, and electron scattering 



To test the SPH-RT method we have examined the collapse 
of a 1 M© molecular cloud of initially uniform density and uni- 
form temperature. The collapse proceeds almost isothermally 
until the density in the centre rises above p ~ 10" 13 gem -3 . 
Then the temperature rises rapidly, the thermal pressure decel- 
erates the collapse, and the first core is formed. The first core 
grows in mass, and contracts quasistatically. When its temper- 
ature reaches T ~ 2, 000 K, the H 2 starts dissociating and the 
second collapse starts, resulting in the formation of the second 
core, i.e. the protostar. Our method reproduces well the results 
of the detailed simulation of Masunaga & Inutsuka (2000): the 
first and the second cores form at similar densities, having sim- 
ilar sizes, and at similar times after the start of the collapse. 
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Fig. 21. Dispersion relation for the thermal relaxation mode. 
The relaxation rates (in units of y) of all SPH particles at seven 
different instants are plotted against the ratio K Q /k (dots that 
saturate to form lines). The theoretical values (calculated from 
Eq.[49]) are also plotted (red line and squares). 

We have also performed the Boss & Myhill (1992) and 
Boss & Bodenheimer (1979) tests, and obtained results very 
similar to those of Whitehouse & Bate (2006). Finally, we 
have performed the thermal relaxation test of Masunaga et al. 
(1998). The geometries treated in this paper establish the fi- 
delity of the method in treating both spherical and flattened 
geometries. Furthermore, the method also performs well on the 
Hubeny (1990) test, which deals with equilibrium discs, and 
hence it can also be applied to disc simulations. We will dis- 
cuss the Hubeny (1990) test, and applications of this method 
to discs in a forthcoming paper (Stamatellos & Whitworth, in 
preparation). 

The new SPH-RT method performs very well, and most im- 
portantly it is very efficient. The computational time is almost 
the same as (only ~ 3 % longer than) an SPH simulation us- 
ing a barotropic equation of state. The method is inherently 
three dimensional, and so it can be used to treat a variety of as- 
trophysical systems, where the radiative processes and thermal 
inertia effects are important. We will report on applications of 
the method in future publications. 
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